# Replication package for 
# "The Economic Leverage of International Organizations in Interstate Disputes"
# Johannes Karreth
# June 30, 2017
# jkarreth@ursinus.edu

# This file: 50_hligo_geo.R
# Purpose: Show spatial and regional patterns of state memberships in IGOs with high leverage
 
rm(list = ls())

# setwd("...")

library("rio")
library("countrycode")
library("ggplot2")
library("maptools")
library("rgdal")
library("countrycode")
library("foreign")
library("RColorBrewer")
library("dplyr")

# Source in functions
source("Functions/theme_jk.R")

# Data
igo_num <- import("igo_overlap_num.csv")
igo_names <- import("Sources/IGO_names_v2.3.csv")
igo_lev3_count_use <- igo_names[igo_names$ionum %in% igo_num$igo_lev3_count_use, ]$IONAME

hligos <- igo_num[is.na(igo_num$igo_lev3_count_use) == FALSE, 1]

members <- rio::import("Sources/IGO_stateunit_v2.3.csv")
hligo_members <- members[ , names(members) %in% igo_lev3_count_use]
hligo_idvars <- members[, c("ccode", "country", "year")]
hligo_members <- cbind(hligo_idvars, hligo_members)

hligo2000 <- hligo_members[hligo_members$year == 2000, -10]
hligo1992 <- hligo_members[hligo_members$year == 1992, c("ccode", "EEC")]
hligo_sum <- merge(x = hligo2000, y = hligo1992, by = "ccode", all.x = TRUE, all.y = FALSE)

hligo_sum$continent <- countrycode(sourcevar = hligo_sum$ccode, 
                                                origin = "cown",
                                                destination = "continent")
hligo_sum[hligo_sum$ccode == 345, ]$continent <- "Europe"
hligo_sum[hligo_sum$ccode == 713, ]$continent <- "Asia"
hligo_sum$memberships <- rowSums(hligo_sum[, 4:20], na.rm = TRUE)

# Map of countries' HLIGO memberships

# This chunk must be executed by itself before the rest of this script
reg <- map_data("world")
P1 <- ggplot(reg, aes(long, lat)) + 
  geom_polygon(aes(group = group), show.legend = FALSE)
worldmap <- map("world", fill = TRUE, col = "transparent", plot = FALSE)
IDs <- sapply(strsplit(worldmap$names, ":"), function(x) x[1])
world_sp <- map2SpatialPolygons(worldmap, IDs = IDs, proj4string=CRS("+proj=longlat +datum=WGS84"))

gpclibPermit()

# Load shapefile

# Might need to set the WD to the Sources folder
# setwd(".../Sources")

world.map <- readOGR(dsn = "TM_WORLD_BORDERS_SIMPL-0.3", 
                     layer = "TM_WORLD_BORDERS_SIMPL-0.3")

# setwd("...")

# Remove Antarctica
world.map <- world.map[world.map$NAME != "Antarctica", ]

# Prepare for ggplot2
world.ggmap <- fortify(world.map, region = "NAME")

# Just to be sure, repeat
world.ggmap <- world.ggmap[world.ggmap$id != "Antarctica", ]

# Supply COW codes
world.ggmap$cowcode <- countrycode(world.ggmap$id, origin = "country.name", destination = "cown")

# See if any COW codes are missing
sort(unique(world.ggmap[is.na(world.ggmap$cowcode) == TRUE, ]$id) )
world.ggmap[world.ggmap$id == "Korea, Democratic People's Republic of", ]$cowcode <- 731
world.ggmap[world.ggmap$id == "Serbia", ]$cowcode <- 345

# Read in HSIGO data by country (prepared with 3 vars for 1950, 1975, 2000)
hligo_sum$cowcode <- hligo_sum$ccode

# Use plyr::join to merge
map.data <- plyr::join(world.ggmap, hligo_sum, by = "cowcode")

# Data frame for plotting
df <- data.frame(id = unique(map.data$id), 
                 cowcode = tapply(map.data$cowcode, map.data$id, mean), 
                 memberships = tapply(map.data$memberships, map.data$id, mean))

# Delete countries w/out COW codes
df <- df[is.na(df$cowcode) == FALSE, ]

# Alternative would be to use scale_color_gradient(low = "white", high = "blue", na.value = "gray")

# Make the plot for mean IGO count 1992+2000
membership.factor <- cut(df$memberships, breaks = c(0, 0.9, 1, 2, 3, 4, 5, 6, 7, 12), include.lowest = TRUE)
p <- ggplot() + geom_map(aes(fill = membership.factor, map_id = id),data = df, map = world.ggmap, color = "black", size = 0.025) + expand_limits(x = world.ggmap$long, y = world.ggmap$lat) + scale_fill_brewer(palette = "Blues") + coord_equal() + theme_classic() + xlab("") + ylab("") +  theme(legend.position = "none", line = element_blank(), text = element_blank(), title = element_blank())

ggsave(p, file = "Output_Tables-and-Figures/hligo_map.pdf", width = 15, height = 10)
  
# Table of HLIGO memberships across continents

hligo_c <- summarize(group_by(hligo_sum, continent),
                     AfDB = sum(AfDB),
                     AsDB = sum(AsDB),
                     CARICOM = sum(CARICOM),
                     ComSec = sum(ComSec),
                     EBRD = sum(EBRD),
                     ECOWAS = sum(ECOWAS),
                     EEC = sum(EEC, na.rm = TRUE),
                     EIB = sum(EIB),
                     EU = sum(EU),
                     IBRD = sum(IBRD),
                     ICfO = sum(ICfO),
                     IFAD = sum(IFAD),
                     IMF = sum(IMF),
                     MIGA = sum(MIGA),
                     Mercosur = sum(Mercosur),
                     SADC = sum(SADC),
                     UEMOA = sum(UEMOA))

t(hligo_c)
hligo_c$represented <- apply(hligo_c[, 2:17], 1, function(x) sum(ifelse(x > 0, 1, 0)))
hligo_c$median <- apply(hligo_c[, 2:17], 1, function(x) median(x))

# Table of HLIGO memberships across regions

hligo_sum$region <- countrycode(sourcevar = hligo_sum$ccode, 
                                             origin = "cown",
                                             destination = "region")
hligo_sum[hligo_sum$ccode == 345, ]$region <- "Southern Europe"
hligo_sum[hligo_sum$ccode == 713, ]$region <- "Eastern Asia"

hligo_r <- summarize(group_by(hligo_sum, region),
                     AfDB = sum(AfDB),
                     AsDB = sum(AsDB),
                     CARICOM = sum(CARICOM),
                     ComSec = sum(ComSec),
                     EBRD = sum(EBRD),
                     ECOWAS = sum(ECOWAS),
                     EEC = sum(EEC, na.rm = TRUE),
                     EIB = sum(EIB),
                     EU = sum(EU),
                     IBRD = sum(IBRD),
                     ICfO = sum(ICfO),
                     IFAD = sum(IFAD),
                     IMF = sum(IMF),
                     MIGA = sum(MIGA),
                     Mercosur = sum(Mercosur),
                     SADC = sum(SADC),
                     UEMOA = sum(UEMOA))

t(hligo_r)
hligo_r$represented <- apply(hligo_r[, 2:17], 1, function(x) sum(ifelse(x > 0, 1, 0)))
hligo_r$median <- apply(hligo_r[, 2:17], 1, function(x) median(x))

hligo_r$regionf <- factor(hligo_r$region, levels = c("Southern Europe", "Eastern Europe", "Western Europe", "Northern Europe", 
                                                     "Melanesia", "Polynesia", "Micronesia", "Australia and New Zealand", 
                                                     "South-Eastern Asia", "Southern Asia", "Eastern Asia", "Central Asia", 
                                                     "Western Asia", "South America", "Central America", "Caribbean", 
                                                     "Northern America", "Southern Africa", "Western Africa", "Middle Africa", 
                                                     "Eastern Africa", "Northern Africa"))
hligo_r$continent <- c("Asia", "Americas", "Americas", 
                     "Asia", "Africa", "Asia", "Europe", 
                     "Asia", "Asia", "Africa", "Africa", 
                     "Americas", "Europe", "Asia", "Americas", 
                     "Asia", "Africa", "Asia", "Europe", 
                     "Africa", "Asia", "Europe")

p <- ggplot(data = hligo_r, aes(x = represented, y = regionf, color = continent)) + geom_point() + geom_segment(aes(x = 0, xend = represented, yend = regionf)) + scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 10)) + theme_jk() + xlab("Number of IGOs with high leverage present in region") + ylab("") + guides(color = FALSE)

ggsave(p, file = "Output_Tables-and-Figures/hligo_region.pdf", width = 6, height = 4)